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Abstract: In this paper we investigate the benefits of using Z(2) ® ^(2) single timesHce 
stochastic sources for the calculation of hght quark physics on the lattice. Meson 2-point 
correlators measured using sources stochastic in only spin and those stochastic in both spin 
and colour indices are compared to point source correlators on the unit gauge and on a 
16^ X 32 Domain Wall QCD ensemble. It is found that the use of stochastic sources gives a 
considerable improvement in statistics for the same computational cost. The neutral kaon 
mixing matrix element Bk is also calculated on this ensemble with stochastic sources, but 
we conclude that the stochastic method offers no significant advantage over the traditional 
gauge-fixed wall source approach which already offers an exact volume average. We also 
discuss the application to semileptonic form factors in conjunction with partially twisted 
boundary conditions. 
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1. Introduction 

Meson correlation functions are fundamental to phenomenological applications of lattice 
QCD. Pseudoscalar states in particular are used to determine quark masses, the LECs of 
the Chiral Effective Lagrangian, and CKM relevant observables such as Bk and the Ki^ 
form factor. 

In this paper we discuss two and three-point functions using the interpolating operator 

Oi,2 = ^Plr^P2 (1.1) 

and its conjugate to create and annihilate mesonic states. Here F is a product of gamma 
matrices set to give the operator the correct quantum numbers, and ipi are quark fields of 
flavour i. 

After Wick contraction, the correlation functions contain a product of spin matrices 
with quark propagators Q. The propagators obey 7^-hermiticity 

g,ix,T ^ y,t) = j^gl{x,T ^ y,th'' , (1.2) 
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for which the arrow indicates the direction of quark flow, and it is understood that the 
source indices of the conjugate propagator lie on the left and those of the unconjugated 
propagator on the right. Even on a lattice of modest size, the propagator matrix is ex- 
tremely large and thus only a subset can be calculated within a reasonable timescale. This 
is performed by solving the matrix equation 

i;{y,t) = Y,M'\y,t ; f,r)r?(f,r) (1.3) 

X,T 

= Y.giy,t^x,T)r]{x,T) (1.4) 

for the 'solution' vector ip, where is a complex vector 'source' occupying some region of 
space, M is the Dirac matrix, and the matrices are contracted over spin and colour indices. 
This equation can be solved using, for example, the iterative conjugate gradient algorithm. 

Typically the point source is used, consisting of unit spin and colour vectors on a single 
space-time point (xq, to). The 12 possible spin and colour source vectors are usually written 
as a unit spin and colour matrix, forming a matrix source fj 



fj{x,t) = I4x4 <8)l3x3 

= 



{x,t) = {xo,to) 
otherwise ' 



where InxN is the N x N unit matrix. 

Solutions evaluated from these sources are matrices consisting of the subset of elements 
of the propagator from a single space-time point to all other points on the lattice, for all 
combinations of spin and colour indices at source and sink, thus requiring 12 inversions of 
the Dirac matrix. These solutions are typically referred to as one-to-all propagators. 

Use of a localised source increases the sensitivity of the measurements to local fluc- 
tuations in the gauge fields; for example propagators evaluated from regions in which the 
Dirac matrix has a localised near zero mode will produce large outliers in the results. The 
statistical distribution of observables should be much better behaved if a volume average 
is included in each measurement. 

This paper is concerned with stochastic vector sources, for which the elements of the 
source are randomly drawn from a distribution T> that is symmetric about zero. A set of 
iVhits randomly generated lattice volume filling sources 

{v^''\x)aaeV\n = l...N^n.}, (1.6) 

referred to as a set of 'hits' of the stochastic source, for which a is a colour index and a a 
spin index, has the property that, in the limit of A^hits ~^ oo 

n=l 

Lattice volume stochastic sources have been used in the past in order to estimate the entire 
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propagator matrix |3|, ^, |6|, |9| . Here the solutions are referred to as stochastic all-to-ah 
propagators. Dong and Liu demonstrated that sources with Z(2) noise V = Z(2) = 
{+1, —1} or generally T) = 'L{N) for any A^, deviate less from the orthonormality condition 



of equation (1.7) for a fixed number of hits than those estimated with Gaussian or 'double- 
hump' Gaussian- like distributions. Foster and Michael H] suggest that the optimal choice 
is the c- number distribution V = Z(2) ® Z(2), which contains random Z(2) numbers in 
both its real and imaginary parts, i.e. 

^ = {^(±l±.)}. (1.8) 

Stochastic all-to-all propagators are generally very noisy, and thus are often abandoned 
in favour of the traditional one-to-all propagators apart from situations in which they are 
necessary (P, ^, such as when the number of gauge configurations is limited and one must 
extract as much information as possible from each. 

This paper details an exploration into an alternate use of stochastic sources for the 
calculation of meson correlators at zero momentum, based upon the work of Foster and 
Michael |^] (appendix), a form of which is referred to as the 'one-end trick' This method 
has been used by the ETM collaboration for the calculation of meson two-point |10, 13, 18 1 



and three-point functions |16|. The aim of this paper is to determine whether these meson 
correlators can be calculated more cheaply and with better statistics than the traditional 
point source, and to investigate the competitiveness of extensions of this method to a range 
of matrix elements compared to the respective traditional approaches. 

The layout of the paper is as follows. The method of the one-end trick is introduced 
in the context of the two-point correlation functions, based upon Foster and Michael's 
|5[ description, followed by the details of the two stochastic source types we have chosen 
to overcome the highlighted issues. The stochastic two-point correlators are compared to 
point source correlators on the unit gauge in order to prove the correct behaviour of the 
stochastic correlators on a trivial gauge. Results for pseudoscalar and vector correlators, 
analysed on a 16'^ X 32 RBC/UKQCD Domain Wall QCD ensemble, are then presented 
and discussed. We then assess the performance of the stochastic source technique in the 
computation of the neutral kaon bag parameter Bk- To this end we compare the evaluation 
of the relevant three-point functions with a single stochastic wall source fixed to Goulomb 
gauge and a naive gauge fixed wall source. We also compare the single-wall approach to the 
two- wall approach of Antonio et al. |17|. We finish with a discussion of the calculation of 
hadronic form factors using the stochastic method, as adopted by the ETM collaboration 



|]J and the RBC & UKQCD collaboration |19|. 



2. Two-point correlation functions 

Using the interpolating operator Oi^2 of equation ( |1.1[ ), the two-point meson correlator is 
defined as 

C{t'-p) = J2e-^''-^y-^H02Mt)OiA^,r)) , (2.1) 

x,y 
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where t' = t — t and p is a momentum. Performing the Wick contraction and neglecting 
the disconnected contribution that appears when ■0i = 'i/'2j we find 

Cit';p) = ^ e-^P<y-^^ X tr {j'Tg,{y, t ^ x, T)Tj'gl{x, r ^ y, t)) , (2.2) 

where the trace is over spin and colour indices and Qi is a quark propagator of flavour i. 

2.1 The One-End Trick 

Consider the meson two-point correlator of equation (|2.2| ) at zero momentum and insert a 
product of Kronecker delta functions 



i^^')xpOl^^^,Jz,T^y,t). (2.3) 



x,y,z 



Here Greek letters represent spin indices and Roman letters colour indices. Using stochastic 
'wall' sources 

t = T 
t^T 



(2.4) 



we can replace the delta functions by a hit average 



SKxScdSx,z = {v''c^\x,T\T)r]l^^\z,T\T))n , 



(2.5) 



using the orthonormality condition of equation (1.7). Inserting this into equation ( p.3p we 
find that the correlator becomes the scalar-product of two solution vectors, one of which 
is dependent on the matrix T: 



(2.6) 



Here 



V^i")(y,t|r) ^Y.g,{y,t^ x,t) r/W(f,^|^) , 

X 



(2.7) 
(2.8) 



and the vectors 7^r^/> and ^ are contracted at the sink location. The source indices are 
contracted automatically by the stochastic average, completing the trace in the A^hits — * oo 
limit. 

The advantages of the above method for the calculation of the entire meson spectrum 
are reduced by the necessity to calculate for each of the 16 T matrices, requiring 16 
inversions per stochastic hit. This can be reduced to 4 inversions per hit by calculating the 
spin structure explicitly. We refer to these as spin-explicit or SEM sources following ref. 
01. These sources are further discussed in section 
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2.2 Pseudoscalar Z2PSWall source 



For the pseudoscalar case F = 7^, the solution vectors ipi ^'^^ V'f ^'^^ identical, allowing for 
calculation of the pseudoscalar meson correlator with only a single inversion of the Dirac 
matrix per hit and valence mass. This completely stochastic source type, defined as per 
equation ( |2.4| ), will henceforth be referred to as a Z2PSWall source. 

Let us consider the structure of the Z2PSWall source in more detail. The spin and 
colour space components can be represented as a single 12-component column vector 
such that the source has the form 

r/(")(x,t|r) ={(^t,4®H(")(x). (2.9) 

The elements of H^") (x) are stochastically sampled from the chosen distribution 

{r!f\x)eV\i = l...l2] (2.10) 

and thus the vectors obey the orthonormality condition 

M(f,y) ^ (HW(x)0SWt(y)>^ ^<5,^^^^Ii2,i2, (2.11) 

where ® is the vector direct product 

M,,-(x,y) ^ (Hf)(f)H;(")(y)>„ = 1...12. (2.12) 

From this we see that M is Hermitian 

m}. (x, y) = M*, (y, x) = M^j (x, y) . (2.13) 

By equation (2.12) we see that the diagonal elements of M (i.e. those with both x = y and 
i = j) are unity. Subtracting these elements, we define the stochastic noise matrix K as 

K{x, z) = Mix, z) - 5^,,- Ii2xi2 . (2.14) 



Thus, we can split the correlator of equation (2^) into signal and noise components 

c{t') ^Y.y{^t\yAr)i^T\yAr))n 

= Cs{t') + AC{t') 



(2.15) 



where 

Cs{t') = Y,tT (gi{y,t ^ x,T)gl{x,T ^ y,t)^ (2.16) 

x,y 

is the gauge-invariant signal component and 

^C{t') = Y,tT (gi{y,t^ x,t)K{x,z)GI{z,t ^ y,t)) (2.17) 

x,y,z 
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is the noise component. Here the trace over sink indices has been reintroduced as this is 
now a product over matrices. The noise component contains a mixture of gauge-invariant 
and gauge-dependent pieces for finite A'hits as well as contributions from other meson 
correlators. This can be seen by decomposing K{x,y) £ C^^^ onto the basis {A^ Tj}, 
composed of the 8 Gell-Mann matrices {A^ I?" = 1 ... 8}, the 3x3 unit matrix Aq = IsxS; 
the 4x4 unit matrix Tq = 1^x4, and the 15 tensor combinations of the gamma matrices 
{Ti \i = 1 . . . 15}. The components of this basis are orthogonal under the trace operation 

tr {\r Fj As Tj) = ttrSrsSij , (2-18) 

where ao = 12 , = 8 |r 7^ 0. Under this decomposition 

Kix,y) = ^Aiix,z) Xr<E)Ti, (2.19) 

i,r 

with c-number coefficients A\.{x,z). Applying this decomposition to the noise component 
AC(t'), we find 

^C{t') = ^^A^^{x,z)iT(gi{y,t^x,T)\r(^Tigl{z,T^y,t)^ . (2.20) 

With reference to equation (^]^), we see that for all components bar the unit matrix 
contribution Aq ® Fq, the spin and colour matrices at the source location are different from 
those at the sink. Therefore these components are formed from the Green's function of 
two different interpolating operators: the pseudoscalar 0{\q 7^) at the sink with the 
polluting 'unwanted operators' 0{\r '^i'y^) at the source. The contaminating noise is 
small in the case of the scalar, vector and tensor Dirac structures as the pseudoscalar is the 
lightest state. These contributions are eliminated in the ensemble average due to parity, 
and also in the A^hits — > 00 limit. The overlap with the axial state Aq is eliminated in the 
hit limit and is empirically smaller in magnitude than the pseudoscalar signal. We shall 
introduce spin-explicit sources for use with non-pseudoscalar measurements. The effects of 
the gauge-dependent terms with x ^ z, which we refer to as 'cross-terms', and components 
with Xr 7^ I3x3 are discussed further in section ^7^ . 



The Z2PSWall source can be implemented within a software framework designed for 
12 X 12 matrix sources, such as the point source discussed in section ||, allowing for the 
reuse of existing propagator contraction code without further modification. This can be 
achieved by placing the stochastic source vector ^("^(x) on the first column of an empty 
12 X 12 matrix <I>(x, t) on each lattice site of the wall 



^>(")(x, t) 



/ : 
H(")(x) 
V : : : : 



t = T 



(2.21) 
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The solution ip' ^'^\y, t\T) is matrix valued 

V^'(")(y,t|r)AC^E-^^B(y>i; ^>^)4c(^>^)> (2-22) 

but with all columns zero bar the first (C = 0). Here A,B,C are spin-colour indices. 
Our inverter, of course, checks for a zero norm source vector before inversion. With this 
implementation, the direct product that forms the stochastic matrix M of equation ( 2.11] ) 
simplifies to the stochastic average of the matrix product 

M{x,y) = mx)^Hy))n, (2.23) 



such that the meson two-point function of equation ( |2.6| ) becomes simply a trace over a 
product of matrices 



C{t';0) = 5^tr (4")(y,t|r)4") ^(y,t|r))„ , (2.24) 



y 

which has the same form as the standard point source meson correlator contraction. 



2.3 Spin-explicit Z2SEMWall sources 

As the Z2PSWall source can only be used for pseudoscalar correlators, we require a source 
type of more general use. We can stochastically estimate the general meson two-point 
correlator of equation (^^) with four inversions by calculating the spin structure explicitly, 
using stochastic noise in colour space only. These sources were used by Viehoff et al. Q for 
the calculation of the matrix element of the axial vector current between proton states, and 
later by Boucaud et al. |1C, ^] for meson correlation functions, where they are referred to 



as 'linked sources'. 

Similarly to the point source, the four spin vectors can be combined into a single 4x4 
unit spin matrix. A different stochastic colour vector ^^"^ is used on every site of the 
timeslice allowing us to retain the spatial and colour delta functions in the hit average, 
while an explicit Kronecker delta is used for the spin components. The source has the 
structure 

77(")(f, t\T) = {5t,r} ® I4X4 ® &\x) , (2.25) 

where 

{it\x)(^V\a = l...?,} (2.26) 
obeys the orthonormality condition 

M(f,y) ^ (e(")(x) ® e^") t(y)>^ ^ 5,^,^^13,3 . (2.27) 

For a finite number of hits, the matrix M can again be decomposed onto the basis of Gell- 
Mann matrices and the unit matrix Aq- As before we expect that the gauge dependent 
terms will be suppressed by the ensemble average. 
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As with the Z2PSWall source, the Z2SEMWall can be placed within a 12 x 12 matrix, 
allowing for the reuse of existing measurement code: 



/l o\ 



/ : 

O 

: 



) 



^>(f, t) 



10 



(2.28) 



10 



\0 1/ 



at t = T and zero elsewhere. 

2.4 Hit averaging and the ensemble average 

The elements of the matrix K are either stochastically generated or zero. Therefore, 
provided the distribution T> is symmetric about zero, the combined probability distribution 
of the gauge configurations U and the stochastic matrix K has the property 



Here we consider the noise components of the correlator containing gauge-dependent terms. 
Terms that break gauge invariance are naturally suppressed by averaging over gauge- 
equivalent configurations within an ensemble. However the Monte Carlo sampling of a 
gauge orbit is typically quite slow as one increases the ensemble size, due to autocorre- 
lations between configurations. The stochastic method improves upon this by explicitly 
removing these terms through the symmetric fluctuations of K about zero. 

If sufficiently many hits per configuration are sampled then the cancellation will be 
near exact, while for a smaller number of hits this will take place stochastically as the 
distribution of gauge fields and sources is jointly sampled. For a large enough ensemble 
the difference between having few hits and having many is likely to be small. 

2.5 Demonstration on the unit gauge 

Following Foster and Michael Q, we use V = Z(2) ®'L{2) noise. We perform our initial 
analysis on the unit gauge configuration, for which all gauge links are unity, in the Domain 
Wall QCD framework. Due to the translational invariance of this gauge field, the point 
source solution is exactly equal to the volume averaged propagator, and therefore this 
configuration is ideal for the demonstration of the convergence of the stochastic correlators 
to the exact volume averaged solution as the number of hits is increased. 

In order to demonstrate the convergence as a function of the number of hits, we 
calculated the pseudoscalar meson two-point correlator with a valence quark mass of 0.04 
in lattice units using different numbers of hits of both Z2PSWall and Z2SEMWall sources. 
For each choice of the number of hits A^hits! we repeat the process 80 times and estimate 
the standard error on the mean of these 80 correlation functions, each comprising of A^its 
solutions. In order to avoid having to generate new data for each A^its) the stochastic 
estimates used for each of the 80 measurements were randomly drawn from a large pool of 
single-hit stochastic estimates of the correlation function, ensuring that for a given Ahits no 




(2.29) 



8 




12 16 20 24 28 32 36 40 44 48 



Figure 1: Demonstration of the dependence of tlie trivial gauge pseudoscalar meson correlator 
on the number of stochastic hits A'hits of Z2PSWall and Z2SEMWall. These are compared to the 
point source correlator which is the exact solution for this gauge configuration. We use 80 separate 
measurements, each consisting of the average of iVhits stochastic estimates in order to determine 
the error. For each iVhits the stochastic estimates are randomly drawn from the pool of estimates 
in order to minimise correlations between the data points. We ensure that each data point is not 
drawn more than once per A'hits such that the measurements are independent. 



data point is drawn more than once. This provides us with 80 independent measurements 
for each N\ii^^ and minimises correlation between values of A'^tiits. 

Figure |^ shows a plot of the means and standard errors of these distributions using 
the correlator at t = 16. For the pseudoscalar two-point function it appears that the 
Z2PSWall is converging better, especially given that it requires one quarter of the number 



of inversions needed for the Z2SEMWall source. However, as described in section 2.3, we 



find that the Z2SEMWall benefits from its exact spin structure in the full calculation. 
2.6 Two-point meson correlator results 
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Lattice size 


16^ X 32 


Action 


Domain Wall 


Gauge Action 


Iwasaki 


Domain wall height 


1.8 


Ls 


16 




2.13 


I (GeV) 


1.73(3) 


Sea quark masses (latt. units) 


ruu = 0.01, rus = 0.04 



Table 1: Ensemble properties. 
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Figure 2: Pseudoscalar meson effective mass plot from averaged correlators with a bin size of • 
configurations. This is not a cost comparison. The points have been slightly shifted for clarity. 
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In order to take the study further, we formed the pseudoscalar meson two-point func- 
tion on 392 configurations (separated by 5 molecular dynamics time-steps) of an RBC- 
UKQCD 16^ X 32 2+1 flavour Domain Wall QCD ensemble with properties detailed in 
table |l|, for which previous results are available for comparison . On each configuration 
the propagators were calculated from Z2PSWall sources with one hit on 12 different source 
timeplanes. For Tgrc = we included 3 further hits. In addition we generated Z2SEMWall 
propagators from 4 timeplanes and point source propagators from two source origins. The 
propagators were calculated using the conjugate gradient algorithm with a residual of 10~^ 
and a valence quark mass of 0.04 in lattice units. This large valence quark mass was chosen 
as it is cheaper to invert; thus allowing for better statistics for a given computational cost. 

In order to take account of autocorrelations in molecular dynamics time, we binned 
over adjacent configurations. Based upon our analysis (appendix ^ and the previous 
analysis |11] we chose a bin size of 40 molecular dynamics time units (8 configurations). 

Figure ^ is an effective mass plot of the averaged correlators with a bin size of 8 
configurations. As before we exclude three of the four Z2PSWall hits on timeslice zero. 
From this figure we chose a constant fit range of 10 — 16 for point source correlators and 
11 - 16 for Z2PSWall and Z2SEMWall correlators. 

Table § contains the results for the pseudoscalar meson mass fits for the various sources 
over all 392 available configurations, where the correlators have been averaged about the 
central timeslice (folded) for better statistics, using the forwards-backwards symmetry of 
the correlator. For some choices of origin or Tgrc the correlation functions show deviations 
from the expected time-dependence which manifests itself in a large value for x^/d.o.f. 
These effects however disappear after averaging the correlation functions over source posi- 
tions. 

The Z2PSWall fitted masses appear to be consistently lower than those of the point 
sources, differing by 5a between the 12-source- aver aged Z2PSWall correlator and the point 
source average. This discrepancy is likely to be caused by statistical fluctuations in the 
point source correlators: The 12-source-averaged Z2PSWall pseudoscalar meson result of 
0.4372(9) is in much better agreement than the 2 point source averaged value of 0.4418(12) 
with the mass obtained in the previous analysis of 0.438(3) [^l|. The central value and 
error estimates of this previous result were obtained by averaging over the pseudoscalar- 
pseudoscalar and axial-axial correlators for several point source smearings and locations 
and the error was scaled by a factor of 1.5 to account for fluctuations in the gauge flelds, 
and as such should not be compared unfavourably with our result. 
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Figure 3: Pseudoscalar effective mass plots at a fixed cost of 4704 inversions of the Dirac matrix. 
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Table 2: Pseudoscalar meson mass fits for the various sources, fitting to range 10 — 16 (point) or 
11 — 16 (stoch. sources), with a bin size of 8 configurations over an ensemble of 392 configurations. 
Nsrc is the total number of sources used in the fit, with the equivalent cost in inversions of the Dirac 
matrix detailed in the next column. The third column of the stochastic source tables contains the 
number of source timeslices used Nt-; the fourth a list of these times Tstc\ and the fifth the number 
of hits (stochastic samples) on those timeslices (A^hits/-^r)- The four independent hits of Z2PSWall 
are distinguished by a Roman letter. 
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At a fixed cost of 12 inversions per configuration (4704 inversions in total), it is evident 
that tlie 12-source averaged Z2PSWall result shows at least a factor of 2 improvement in 
the statistical error over the point source. The 3-source averaged Z2SEMWall results agree 
with the Z2PSWall result, and also show a consistent factor of 2 improvement in errors 
over the point source at the same cost. 

The third- and second-to-last lines in the Z2PSWall section of table || allow for a fixed 
cost comparison between the use of four stochastic hits upon a single source timeslice and 
a single hit on four different timeslices. The 40% reduction in the statistical error suggests 
that one should preferentially choose a new three-volume sample of the gauge field when 
forming a new stochastic source, separated in space-time and molecular dynamics time in 
order to maximise decorrelation between the samples. 

Figure ^ plots the effective mass as a function of time at the cost of 4704 inversions. It 
is evident that both the Z2PSWall and Z2SEMWall source types give significantly better 
plateaus than a single point source. The plateaus for the stochastic sources appear to be 
very similar. We believe this displays a spectacular improvement for pseudoscalar masses 
at no additional cost. Based upon these data we conclude that the difference in the quality 
of the stochastic source results at fixed cost is not large enough to warrant using Z2PSWall 
sources over the spin-explicit Z2SEMWall sources which can be used for a larger number 
of measurements. 

In order to estimate the effectiveness of the Z2SEMWall sources for other measure- 
ments, we consider the vector meson correlator with interpolating operator Oi^2 = '4'i^'^ip2- 
Figure |^ shows the vector meson effective mass for the 4 combined Z2SEMWall sources 
and the 2 point sources, where we have averaged over the three spatial gamma matrix 
correlators and folded about the central timeslice for better statistics. Based upon this 
plot we chose a fit range of 9-16 for both source types. This is not a fixed cost comparison. 

Table |3| shows the results of the fits to the vector meson correlator. As before, at fixed 
cost we see a reduction in error by a factor of around 2 over the point source values. Also, 
as with the pseudoscalar two-point function, the plateau of the vector meson effective mass 
(figure ^) for the stochastic source type is noticeably better than that of the point source. 
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Figure 4: Vector meson effective mass plot from averaged Z2SEMWall and point source correlators 
with a bin size of 8 configurations. This is not a fixed cost comparison, but can be used to select 
the fit ranges. 
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Figure 5: Vector meson effective mass plots at a fixed cost of 4704 inversions of the Dirac matrix. 
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Point 


N 

-"src 


Cost 


Orij 


5in(s) 




Mass 


xVd.oI. 


1 


4704 


Xl = 


E (0,0,0,0) 




0.656(+10)(-9) 


0.404 


1 


4704 


X2 = 


= (8,8,8,16) 




0.657(+8)(-9) 


0.172 


2 


9408 


Xl , 






0.657(+7)(-7) 


0.372 


Z2SEMWall 


Nsrc 


Cost 


Nr 






Mass 


xVd.o.f. 


1 


1568 


1 





1 


0.642(+8)(-8) 


0.207 


1 


1568 


1 


4 


1 


0.660(+10)(-ll) 


0.328 


1 


1568 


1 


8 


1 


0.642(+10)(-9) 


0.407 


1 


1568 


1 


12 


1 


0.637(+9)(-9) 


0.538 


3 


4704 


3 


0,4,8 


1 


0.649(+5)(-5) 


0.153 


3 


4704 


3 


0,4,12 


1 


0.647(+5)(-6) 


0.101 


3 


4704 


3 


0,8,12 


1 


0.641(+5)(-5) 


0.377 


3 


4704 


3 


4,8,12 


1 


0.646(+6)(-6) 


0.457 


4 


6272 


4 


0,4,8,12 


1 


0.646(+4)(-5) 


0.179 



Table 3: Vector meson mass fits for the varfous sources, fitting to range 9 — 16 with a bin size of 
8 configurations. Here we have used the conventions estabhshed in table 0. 
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3. Stochastic calculation of the Bp^ three-point function 



The success of the stochastic method for two-point functions motivates us to consider 
extending the technique to include meson three-point functions. In this paper we consider 
the kaon bag parameter Bk, which describes the mixing between the neutral Kq and Kq 
mesons: 

platt_ {K'^\OvV+Aa\K') .... 

^ - |(i^O|Ao|0)(0|^o|^°) ' ^ ^ 

where \K^) is a neutral kaon state, is the zeroth component of the axial current and 
Ovv+AA is the AS = 2 double weak decay operator responsible for the mixing, which has 
the form 

Ovv+AA = (srd) (srd) + (s7'7^d) (sj^l'^d) . (3.2) 
The Green's function {K^ {ti)\Orr{t)\K'^ {t2)) has two Wick contractions 

W\t;T) = Y,tr My,t;T) X tr ^2{y,t;r) (3.3) 



and 



where 



W\t;T) = j;tr($i(y,t;r)cl>2(y,t;r)) , (3.4) 



^i{y,t;T) =J2Gs{y,t ^ u,U)gl{v,ti ^ y,t)j^T . (3.5) 

u,v 

Here, y is a vector on the interaction timeslice t; T is the relevant product of gamma 
matrices T £ {t'^jT^t'^}; the subscript on the propagator refers to the quark flavour and 
the 7^-Hermiticity of the propagator has been used as before. 

RBC & UKQCD typically use propagators calculated from a pair of spatially separated 
gauge fixed walls of Kronecker delta sources (referred to as GFWall sources). These are 
usually calculated with both periodic (p) and antiperiodic (a) boundary conditions, using 
the p + a combination to eliminate unwanted round-the-world contributions to the three- 
point function by doubling the periodicity of the meson's propagation. This method has 
been used to calculate Bk on this 16^ x32 ensemble |17] using GFWall sources on timeslices 
ti = 5 and t2 = 27. 

In this paper we use a variant of this method, using a single GFWall source at time 
r with any antiperiodic signs implemented on all time-directed links Ut{T — l,x). Here we 
may take the p+a combination as a forwards propagating solution and p — a as a backwards 



propagating solution. This method has been used by Aubin et al. |14] for the removal of 



round-the-world pion propagation in the calculation of the pseudoscalar decay constant. 
3.1 Stochastic Bk 



The form of equation (p^ is suitable for calculation using both Z2PSWall and Z2SEMWall 
sources. However the stochastic cancellation to a delta function in space supresses the cross- 
terms by design. This property can be removed by choosing the same set of stochastic 
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numbers on each site of the timeshce for a given hit, such that M{x, y) becomes position 
independent and yet still retains the delta function in spin and colour space in the large hit 
limit. Thus, we introduce two new source types with this property: The Z2PSGFWall and 
the Z2SEMGFWall source which otherwise share the same source structure as the existing 
stochastic types. These sources should be used on gauge fixed configurations. 

3.2 Bk results 

We calculated Bk on the 16'^ x 32 ensemble detailed in table ||, using a valence quark 
mass of 0.04 lattice units for all propagators. Fixing to Coulomb gauge, we used propa- 
gators generated from the Z2PSGFWall and Z2SEMGFWall stochastic source types and 
the standard GFWall source. We also generated propagators from the Z2PSWall and 



Z2SEMWall types without gauge fixing. Following the discussion in section 2.4 we gen- 
erate a single stochastic hit per configuration for all source types bar the Z2PSWall and 
Z2PSGFWall. For the latter types, the limited size of the ensemble forced us to increase 
the number of hits per configuration to four in order to compare with the GFWall results 
at reasonable statistics. For a fair comparison at a given cost, the configurations used for 
the GFWall correlators were spread across the ensemble in order to reduce the effect of 
autocorrelations . 

Upon loading each gauge configuration, we perform a spatial translation by a pred- 
ermined four- vector d. With every new configuration d is incremented by an amount A, 
allowing us to spread the sources throughout the lattice volume without the need to alter 
the location of the timeplane upon which we apply the boundary conditions: The sources 
are always placed at t = on the shifted configuration with the boundary conditions applied 
on the boundary between t = T — 1 and t = 0. We thus intended that the n*'^ configuration 
be shifted dn = di + {n — 1)A, where the periodicity of the lattice is implicit. While this 
rule was mostly followed in contiguous segments, due to restarting the code this rule was 
interrupted at several points in the chain. The actual source origins are widely distributed 
and for the most part follow the above rule, and thus this will not substantially affect the 
conclusions. Our subsequent production running for phenomenological calculations using 
these methods follow the above rule strictly. 

The fit ranges were chosen based upon the plateau range of the PA correlators in the 
denominator of equation and the errors were estimated using the jackknife procedure. 

The data were binned over a minimum of 40 molecular dynamics time units as before. 
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Source Type 


#conf. 


Fit range 


Fit value 


xVd.o.f. 


Scaled fit 


Z2PSWaU 


384 


7-22 


0.6338(75) 


1.882 


0.6338(160) 






9-24 


0.6736(88) 


1.039 


0.6736(188) 






7-25 


0.6374(83) 


1.185 


0.6374(177) 






9-24 


0.6479(68) 


1.162 


0.6479(145) 


Z2PSGFWall 


384 


7-24 


0.6155(79) 


0.721 


0.6155(198) 






9-23 


0.6595(81) 


1.092 


0.6595(203) 






10-23 


0.6492(94) 


0.857 


0.6492(235) 






8-22 


0.6250(85) 


1.108 


0.6250(213) 


Z2SEMWall 


128 


9-24 


0.6752(56) 


0.753 


0.6752(173) 


Z2SEMGFWall 


128 


10-24 


0.6685(43) 


1.874 


0.6685(67) 



Table 4: Results for Bk on the 16^ x 32 ensemble for the various source types calculated at a 
fixed cost of 384 inversions. The number of configurations is given in the second column. Four 
independent hits over the same set of configurations were calculated for the Z2PS types; here we 
quote the results of independent fits to each of these sets in order to demonstrate the fluctuations in 
the correlators resulting from the choice of different random numbers. These data are inconsistent 
and thus we scale the errors by a PDG scale factor with the results given in the last column. PDG 
scale factors are calculated for the Z2SEM types by splitting the available data into two sets and 
performing separate fits as discussed below. 



Source Type 


Fit range 


Fit value 


xVd.o.f. 


Scaled fit 


Z2SEMWall 


9-23 


0.6626(39) 


1.010 


0.6626(121) 




8-23 


0.6815(49) 


1.193 


0.6815(152) 


Z2SEMGFWall 


8-24 


0.6665(48) 


0.801 


0.6665(75) 




7 25 


0.6572(38) 


1.138 


0.6572(59) 


GFWall 


7-25 


0.6590(28) 


1.278 






8-25 


0.6579(24) 


1.081 





Table 5: Bk fits over 2 sets of 192 configurations with a separation of 10 configurations. The two 
sets are staggered by 5 configurations such that there is no overlap, thus approximating 2 hits on 
the same configurations. 



Source Type 


#conf. 


Fit range 


Fit value 


xVd.o.f. 


Scaled fit 


Z2PSWall 


4 X 384 


8-24 


0.6548(51) 


0.261 


0.6548(109) 


Z2PSGFWall 


4 X 384 


8-25 


0.6365(50) 


0.676 


0.6365(125) 


Z2SEMWall 


384 


7-23 


0.6728(30) 


0.755 


0.6728(93) 


Z2SEMGFWall 


384 


7-24 


0.6653(28) 


0.913 


0.6653(43) 


GFWall 


128 


9-23 


0.6554(32) 


2.000 





Table 6: Results for Bk on the 16^ x 32 ensemble for the various source types. These data are 
calculated at a fixed cost of 1536 inversions, where the Z2PS types were evaluated for 4 hits over 
the same 384 configurations. The value quoted in Antonio et al. ||l^ is 0.659(3). 
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Source Configuration 


Fit range 


Fit value 


xVd.o.f. 


1 wall, t = 


7-25 


0.6591(28) 


1.353 


2 walls, h = 5,t2 = 27 


9-23 


0.6634(26) 


0.779 



Table 7: Fits to the Bk three-point function using GFWall sources at fixed cost, comparing the 
single-wall source method to the traditional two-wall source method. 



Table ^ shows the results for the stochastic types calculated at a fixed cost of 384 
inversions. The GFWall results have been omitted from this table due to the lack of 
statistics at this number of inversions (32 configurations). The fits to each of the 4 hits 
of Z2PSWall and Z2PSGFWall show discrepancies in their central values outside of the 
quoted errors, with a combined x^/d.o.f. of 4.569 and 6.274 respectively. No improvement 
in the agreement was found by increasing the bin size. We therefore apply a PDG scale 
factor of x^/d-o.f . to the Z2PS error bars in order to account for this unlikely, high 
X^/d.o.f. The results of this scaling are given in the far right column. 

In order to investigate the appropriate scaling factors for the other source types, we 
consider the fits to 2 sets of 192 configurations with a separation of 10 and a bin size of 
40 MD time units (4 configurations). The sets are staggered by 5 time units such that, 
due to correlations between nearby configurations, this method approximates two hits on 
the same 192 configurations without the need for further computation. The results of this 
analysis are presented in table |5|. From each of these pairs of fits we calculate the combined 
X^/d.o.f. allowing us to estimate the PDG scale factor as above. It is clear that the GFWall 
results agree very well within errors, with a x^/d.o.f. of 0.091 and therefore need no scaling. 
However, the agreement between the two fits for the Z2SEM correlators is poorer, with a 
xVd-o.f. of 9.591 for the Z2SEMWall results and 2.436 for the Z2SEMGFWall. Again we 
scale the Z2SEM results by the PDG scale factor, with the results included in tables Q and 
i 

Combining all available hits of the stochastic source types at a fixed cost of 1536 
inversions allows for their comparison with the GFWall source correlators calculated on 128 
configurations. These results are presented in table |6|, where we have scaled the stochastic 
source results by their appropriate PDG scale factors as before. After rescaling, all of the 
fits presented agree with the value quoted in Antonio et al. |17] of 0.659(3). From these data 



it is evident that the stochastic approach shows no advantage over the traditional method 
for the calculation of the Bk matrix element, although the Z2SEMGFWall correlators, 
which have the same structure as the GFWall correlators in the large A'^hits limit, give 
comparable results at the same cost. 

3.3 Comparison of the two-wall and single-wall approach to Bk 

Using data from ref. we are able to compare the single-wall approach to the calculation 
of Bk with the traditional two-wall method. Table |^ compares the fits to Bk calculated 
with the two methods at fixed cost, using 196 configurations of single- wall data and 98 
configurations of two- wall data. These sets of configurations overlap, and have configuration 
separations of 10 and 20 molecular dynamics time units respectively such that the set 
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Figure 6: A comparison of tlie Bk plateau calculated using the single-wall and two- wall approaches 
at a fixed cost in inversions. The two-wall sources reside on timeslices 5 and 27. 



lengths are similar. Both sets of data were analysed with a bin size of 40 MD time units. 
The gauge fields of the single-wall calculation were shifted between configurations in order 
to reduce the effects of autocorrelations, whereas the two wall sources were fixed at t = 
5 and 27. Figure ^ compares the Bk plateaus of the two methods. Prom these data 
we conclude that the single-wall method gives equivalent results at fixed cost. However, 
the method affords the sampling of more timeslices and configurations than the two-wall 
approach for the same cost. 



4. Stochastic calculation of three-point hadronic form factors 

The TT^lui form factor Ki^ and the pion electromagnetic form factor are phenomeno- 

logically interesting parameters calculated from relatively simple meson three-point func- 
tions and the meson two-point correlators discussed in section |2|. 

Using the notation of Boyle et al. |12], the three-point functions have the form 

Xf,X 

- ^'^^ {Pf{pf)\v,io)\PM)) 



x|0(i^_t)e-^'(*-*')-^/(*/-*) - 0(t-t;)e-^'(^+*>-*)-^/(*-*/)} , (4.1) 

where pseudoscalar {i, f G {vr, K}) initial states Pj and final states Pf, with energies Ei and 
Ef respectively, are created by the interpolating operators Oij = V'i7^V'2 with fermions of 
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the appropriate flavour. is the vector current operator with appropriate quark flavours 
to ahow the transition, and Zj = = ( | 0/(0, 0)| Py ). The source and sink timeplanes 
ti and tf are typicaUy fixed, with a large time separation to remove the round-the-world 
contribution. 

After Wick contraction, the three-point function becomes the trace over contracted 
propagators 

tr (g(0, ti ^ Xf,tf)j^g{xf,tf ^ X, t)j''g{x, t ^ O, ti)j^^ . (4.2) 

These propagators can be determined from a single source point {0, ti) using a standard 
point source for the propagator g{x,t ^ 0, tj) and a sequential propagator 

g(0, t^^pf,tf^x,t) = Y^ 75 {g{x, t ^ Xf,tfh^g{xf,tf ^ 0, t,) e-*P^/"'/)%5 (4.3) 

for the product G{0,ti <— Xf,tf)^^Q{xf,tf <— x,t), including a Fourier transform over xj 
to momentum pf at the sink timeplane tf. The trace is then Fourier transformed over the 
vertex position x with the phase factor e^^^^~'^f^'^ to complete the three-point correlator. 

For zero spatial momentum pi at ti, a stochastic wall source can be used in place 
of the traditional point source, giving an estimate of the spatial volume average. The 
stochastic averaging to Kronecker deltas will occur on the source timeplane, with the 
second leg of the sequential propagator inverted on the stochastic solution vector. We note 
that although our stochastic sources explicitly project to zero source momentum, partially 
twisted boundary conditions |12, @] can be used in conjunction with this method to apply 



a residual momentum pi. 

From equation (^]^) we see that the propagators are contracted at the source without an 
intervening gamma matrix, and thus these three-point functions are suitable for calculation 
with Z2PSWall sources, as well as the more general Z2SEMWall. 



This method has been adopted by the RBC & UKQCD collaboration |1£] for the 
calculation of meson form factors. In the above, the authors compare the point source 
and stochastic approaches at a fixed statistical accuracy, concluding that the stochastic 
approach is vastly superior to the point source, offering similar statistical errors for less 
than ten percent of the computational cost. A similar approach is also used by the ETM 
collaboration jl6[ |. 



5. Summary and conclusions 

In this paper we have detailed our investigations into the use of stochastic wall sources 
for the calculation of meson two-point and three-point functions using the one-end trick 
fs], P]. We emphasise that this one-end method is a different application of stochastic 
sources to the method of approximating the all-to-all propagator for the calculation of 
disconnected correlation functions: The one-end trick uses the properties of the stochastic 
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sources to offer a volume averaging of tlie standard connected correlation function alongside 
an overall reduction in computational cost. 



In section 2.2 we have described in detail the structure of two Z(2) ^(2) stochastic 
wall source types, namely the Z2PSWall and Z2SEMWall, where the former is random in 
spin and colour space and the latter only in colour space, discussing the form of the noise 
introduced into various measurements. 

The viability of these source types for the calculation of meson two-point functions 
on the unit gauge and on a 16^ X 32 Domain Wall QCD ensemble was demonstrated in 
sections and 2.6. We have shown that both stochastic source types give errors on 
the pseudoscalar meson mass that are smaller by a factor of two or more than those of 
the conventional point source approach at the same cost. In addition to the reduced error, 
there is a substantial improvement in the quality of the plateaus (figure ^) inspiring greater 
confidence in the results. 

In principle we believe that wall source techniques offer better sampling of low probabil- 
ity tails of the QCD functional distribution, for example physically significant contributions 
from rare, low eigenvalue modes of the Dirac matrix. Such modes are likely to produce 
outliers when sampled by a point source. The relative improvement of wall sources over 
point sources will likely increase with lattice volume and decreasing quark mass, a feature 
in common with the low mode averaging approach. 

The Z2SEMWall is also shown to be viable for other meson correlators, showing im- 
proved statistical error over the point source for the vector meson mass. Thus we conclude 
that meson spectrum measurements such as masses and decay constants can be calculated 
with improved precision and confidence using the stochastic method. 

In section ^ we have shown that stochastic sources are viable for the calculation of the 
kaon bag parameter, Bk- We have included two extra stochastic source types in the analysis 
that stochastically estimate the Coulomb gauge fixed wall source {GFWall), each treating 
the spin-colour trace differently. These are referred to as Z2PSGFWall and Z2SEMGFWall. 
However, we found that the more complex structure of these three-point functions is less 
well treated by stochastic methods. Multiple measurements using stochastic sources on the 
same configurations showed disagreements in their central values outside of the jackknife 
error bars, forcing us to apply a PDG scale factor of x^/d-o.f . to the error bars of 
the stochastic results. We conclude that for three-point matrix elements of Oyv+AA-, the 
stochastic method offers no corresponding substantial gain over the traditional GFWall 
method. 

We also find that the use of a single GFWall source calculated with periodic and an- 
tiperiodic boundary conditions, from which we calculate the forwards {p+a) and backwards 
{p — a) propagating components, offers comparable cost-effectiveness to the two- wall meth- 
ods, but may allow more time origins or more configurations to be used when measurement 
cost is the limiting factor. 

Finally, in section |^ we discuss a method of stochastically estimating hadronic form 
factors. This method has been adopted by RBC&UKQCD for a calculation of the Ki-^ form 
factor []T9|, with the conclusion that a significant reduction in computational cost can be 
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achieved for the same statistical error using the stochastic method coupled with partially 
twisted boundary conditions. 

The relative difference in gain between the Ki^ form factor and Bk can be explained 
as follows. The standard GFWall method for the calculation of Bk already provides an 
exact three- volume average of the operator insertion point, and thus the only benefit of the 
stochastic wall method is in the reduced cost of the spin-color tracing in the pseudoscalar 
interpolating operators. Our results suggest this is empirically ineffective. 

In contrast the requirement of non-zero momentum for Ki^ results in a comparison 
of a localised source to a three-volume average, and there is much more scope for the 
stochastic volume average to gain. In this calculation, of course, momentum is injected 
using partially twisted boundary conditions, and cost of requiring multiple inversions for 
different momenta must be included. It is certainly possible that, similar to Bk, a gauge 
fixed wall source in combination with partially twisted boundary conditions could result 
in an even greater improvement than the stochastic wall, by similarly providing an non- 
approximate volume average at similar cost. This is something we intend to study. 
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A. Binning analysis of the 16^ x 32 ensemble 

In order to prove any gain displayed in statistical errors is real and not associated with 
the spurious introduction of additional, correlated measurements we must ensure the in- 
dependence of our extra data points. Thus we bin over adjacent configurations. The bin 
size was chosen by averaging over pseudoscalar meson correlators independently for each 
source type {Z2PSWaU, Z2SEMWall and point) and fitting to a fixed range of 10 — 16 
while varying the bin size. We chose to include only a single Z2PSWall hit on timeslice 
zero, as using all four available hits would weight the statistics in favour of this timeslice, 
making it more difficult to estimate the autocorrelations. 



Figure A-1 shows the change in the estimated error on the fit as a function of the bin 
size for the chosen fit range, given as a fraction of the error on the unbinned data. The lack 
of smoothness of the curves in the figure indicates that the statistics are not good enough 
to provide anything but a rough estimate of the required bin size. The previous analysis 
1 11], using a combination of many source smearings and types, placed a lower bound of 20 
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Figure A-1: Fractional change in the estimated error on the fits to the source averaged pseudoscalar 
meson correlators, as a function of bin size. Fits are performed to the chosen fit range 10 — 16. 

molecular dynamics time units on the integrated autocorrelation length of this ensemble, 
corresponding to a separation of 40 molecular dynamics time units (8 configurations) for 
independent measurements. This is in agreement with our analysis. Based upon this work 
and for comparison with the above we chose a final bin size of 8 configurations. 
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